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Abstract. We present a new algorithm for solving ideal relativistic hydrodynamics based on 
Godunov method with an exact solution of Riemann problem for an arbitrary equation of state. 
Standard numerical tests are executed, such as sound wave propagation and the shock tube 
problem. Low numerical viscosity and high precision are attained with proper discretization. 


1. Introduction 

Ultrarelativstic heavy ion collisions are nowadays extensively studied both in theory and 
experiment. The collision leads to the creation of hot and dense matter, quark gluon plasma, that 
soon thermalizes and expands. This expansion has been successfully described by relativistic 
hydrodynamics EE]. Since the initial conditions used in hydrodynamic modeling are often 
complicated, possibly leading to shock waves, sophisticated numerical methods must be used. 
A class of high-resolution shock-capturing methods particularly suited to solve this problem are 
Godunov methods. We present a new numerical scheme for ideal relativistic hydrodynamics 
using the exact solution of Riemann problem for an arbitrary equation of state (EoS). 

2. Hydrodynamical modeling 

The ideal relativistic hydrodynamic equations have the following form: 

d„n.* = 0, (1) 

d„T^ = 0 , 

where n^ is the flow of a conserved charge and is the energy and momentum tensor in 
the non-dissipative case. In nuclear collisions the relevant charge is the baryon number, thus 
ttT = tibu where u^ = 7(1,#) is the flow velocity, 7 being the Lorentz factor. Our model 
covers nuclear collisions at highest energies, where the net baryon density is practically zero. 
Therefore we do not consider the first equation at all and solve only the second equation that 
expresses the conservation of energy and momentum. The energy and momentum tensor has, 
in the ideal case, this explicit form: 

T|J = (e + p) U V-j«r, 


( 2 ) 


where e is the energy density, p is the pressure and g^ v — diag(l, —1, —1, —1) is the Minkowski 
metric. The second equation can be rewritten in a different form, useful for numerical 
implementation, where we employ a vector of conserved variables U and the spatial flow of 
the conserved variables is F(U)\ 

d t U + d x F(U) = 0, (3) 

where 

u= ((e+pb 2 -f>, (e + pbV, (4) 

(t+p)'y 2 v 2 ,{e+ph 2 v 3 ) T , 

F l = ((e + p)"f 2 v\ (e + p)^ 2 v l v 1 + 8 ll p, (5) 

(e + p) r ) 2 v l v 2 + 8 l2 p, (e + p)-f 2 v l v 3 + <5* 3 p) T . 

The presented equations are solved numerically in the latter form using Godunov method with 
an exact Riemann solution at the interface. In this method we consider a piecewise constant 
distribution of variables inside the cells of a numerical grid. To obtain the value of conserved 
variables at the next time-step, a local Riemann problem is solved at each interface of two 
neighbouring cells. Its solution allows us to compute fluxes of conserved variables F(U) at the 
interface M. For a given cell, we then obtain the values of conserved variables at the next 
time-step by averaging the fluxes on the left and right boundary of the cell: 

Ut+M = u t + ^(F l+1/2 - Fi_ 1/2 ), (6) 

where U\ is the value of a variable in the i th cell at time t and F i _ 1 / 2 (F i+1 / 2 ) is the flux at its 
left (right) boundary. 


3. Numerical tests 

The sound wave propagation test consists of simulating a sound wave of one wavelength in the 
numerical grid. We impose the following initial conditions: 


Pinit ( x ) 
Vinit( x ) — 


c . 27 TX 

= Po + opsin-^-, 

dp . 27 TX 

c s (e 0 + Po) Sm A 


( 7 ) 


with parameters po = 10 3 fm -4 , dp — 10 -1 fm -4 . Since the variation of pressure is sufficiently 
small dp po 5 we can consider the linearized analytic solution, a sound wave of velocity c s , and 
compare this solution to the values obtained by numerical computation to evaluate the precision 
of our numerical scheme. The precision is studied with LI norm evaluated after one time-step, 
that corresponds to the period of the wave, with different numbers of cells in the numerical grid 


Ncell 

N C ell 

L(p(N ce u),p s ) = V \p(Xi\Nceii) - p a (Xi) |—- 

P>cell 


( 8 ) 


The dependence of the LI norm on the number of cells is shown in the left panel of Figjl] The 
precision improves with finer discretization, as expected, and is very good for grids larger than 
N ce ii > 500. Despite the fact that we are using ideal hydrodynamics, the numerical computation 
introduces dissipation into our solution. Since quark gluon plasma is expected to have a very 
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Figure 1. Dependence of LI norm (left) and numerical viscosity r] nurn (right) on number of 
cells in the numerical grid. 


low viscosity, the artificial dissipation of the numerical scheme must be kept very low too. We 
have evaluated the numerical viscosity of the scheme r ] num , using LI norm: 


3A , 

'Hnum — "F ^0/ ^ 


1 


7r 


2 A Sp 


L(p(N ce ii),p s ) 


(9) 


We present the dependence of numerical viscosity r] nurn on the number of cells in the grid in the 
right panel of Fig|lJ Similarly to the dependence of LI norm, it decreases with number of cells 
and its values are adequately small. We also evaluated a more suitable parameter - the ratio of 
numerical viscosity and entropy density in our scheme r] nurn /s. We present its values in Figj2] 
together with values of r]/s for pion gas [6] and the limiting value of r]/s for quark-gluon plasma 
from AdS/CFT calculations 1/47r [7]. From this comparison we can see that the ratio r] nurn /s 
is sufficiently small for modelling of quark-gluon plasma and will not influence the results due 
to numeric effects more than the physical viscosity. For simulations in ID this test shows that 
suitable grid sizes are N ce u > 400. The 3D case remains to be checked. 



Figure 2. Dependence of the numerical viscosity to entropy density ratio r] nurn /s (blue points) 
compared to the limiting value r\/s — 1/47r (red) and rj/s of pion gas (green band). 

We continue with the shock tube problem. It consists of imposing special initial conditions 
with a discontinuity in energy density in the middle of the spatial domain between two constant 










states. The initial conditions in energy density are = 16 GeV (e# = 1 GeV) temperature in 
the left (right) half of the numerical grid. The initial velocity is zero over the whole grid. With 
time, we expect to see the dissolution of the discontinuity into a rarefaction wave propagating 
to the left, to the region of higher energy density, and a shock wave propagating to the right, 
where energy density is lower. The shock tube problem has an analytic solution, which allows a 
good comparison between the numerical and exact solution. In the left panel of Figj3]we show 
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Figure 3. Profile of the energy density e (left) and the velocity v (right) in the numerical grid 
after 100 time-steps (our scheme in blue, analytic solution in black). 

the profile of energy density in the grid after 100 time-steps. We see that the numerical solution 
(blue) is comparable to the analytic solution (black). The scheme is able to handle the initial 
discontinuity very well. The profile of velocity after 100 time-steps, which we show in the right 
panel of Figj3j replicates the analytic solution similarly well. 

4. Outlook 

We have built and tested an ideal relativistic hydrodynamic scheme based on the exact solution 
of Riemann problem. The presented tests in one spatial dimension show a good resolution 
and ability to capture shock and rarefaction very well. We will extend this scheme to three 
dimensions and then apply it in description of the flow in ultrarelativistic nuclear collisions. 
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